############################################################################
# simulation function that estimates deforestation and emission
# with a carbon price
pred_fun_defor <- function(mdl_loss, 
                     mdl_gain,
                     dataf, 
                     forest_cv,
                     ag_rev_name, 
                     carbon_name, 
                     carbon_price,
                     biome_cat_df,
                     start_for_cv_bin,
                     npv_carbon_refor_df,
                     plant_per_df,
                     soil_carbon_refor_df,
                     npv_soil_carbon_refor_df){

  
  #mdl_loss = cont_loss_list
  #mdl_gain = cont_gain_list
  #dataf = cont_df
  #forest_cv = start_for_cv
  #ag_rev_name = ag_rev_name 
  #carbon_name = carbon_name 
  #carbon_price = 20
  #start_for_cv_bin = start_for_cv_bin
  #biome_cat_df = biome_cat_df
  #npv_carbon_refor_df = npv_carbon_refor_df
  #plant_per_df = plant_per_df
  #soil_carbon_refor_df = soil_carbon_refor_df
  #npv_soil_carbon_refor_df = npv_soil_carbon_refor_df
  
  
  #### gain prediction 
  # create a variable of character variables 
  bio_fix = paste0("F_biome_num", dataf$biome_num_f)
  bio_fix = mdl_gain[match(bio_fix, names(mdl_gain))]
  
  # create continent variable 
  cont_fix = paste0("F_cont_num_f", (ifelse(dataf$cont_var == "LatinAmerica", 1,
                                            ifelse(dataf$cont_var == "Africa", 2, 3))))
  cont_fix = mdl_gain[match(cont_fix, names(mdl_gain))]
  
  # calculate biome reforestation category 
  biome_npv_cat = biome_cat_df$biome_cat[match(dataf$biome_num, biome_cat_df$biome_num)]
  cont_cat = ifelse(dataf$country == "Brazil", 4,
                    ifelse(dataf$country == "Cambodia", 5,
                           ifelse(dataf$country == "Colombia", 6,
                                  ifelse(dataf$country == "Indonesia", 7,
                                         ifelse(dataf$country == "Liberia", 8,
                                                ifelse(dataf$country == "Malaysia", 9,
                                                       ifelse(dataf$country == "Peru", 10,
                                                              ifelse(dataf$cont_var == "LatinAmerica", 1,
                                                                     ifelse(dataf$cont_var == "Africa", 2, 3)))))))))
  
  match_df = data.frame("uni_ID" = dataf$uni_ID,
                        "biome_npv_cat" = biome_npv_cat,
                        "cont_cat" = cont_cat, 
                        stringsAsFactors = FALSE)
  match_df = left_join(match_df, plant_per_df, by = c("cont_cat"= "contin_cat"))
  
  npv_refor_val = left_join(match_df, npv_carbon_refor_df, by = c("biome_npv_cat" = "biome_category",
                                                                  "cont_cat" = "contin_cat"))
  npv_refor_val$ind_var = ifelse(is.na(npv_refor_val$biome_npv_cat), 0, 1)
  
  
  # create soil carbon reforestation data 
  soil_df = data.frame("uni_ID" = dataf$uni_ID,
                       "cont_cat" = cont_cat, 
                       stringsAsFactors = FALSE)
  
  soil_df = left_join(soil_df, npv_soil_carbon_refor_df, by = c("cont_cat" = "match"))
  
  # carbon gain 
  carbon_gain = ((npv_refor_val$discount_rate_10per * 1.26 * 3.67) + 
                   (dataf$soil_crb_ref * soil_df$discount_rate_10per)) * (carbon_price / 9.567)
  
  # calculate prediction model 
  pred_gain = mdl_gain["(Intercept)"] + 
    mdl_gain[ag_rev_name] * (dataf[,ag_rev_name] - carbon_gain)  +
    dataf$elev * mdl_gain["elev"] + 
    dataf$slope * mdl_gain["slope"] + 
    dataf$iucnlow_2000 * mdl_gain["iucnlow_2000"] + 
    dataf$iucnhigh_2000 * mdl_gain["iucnhigh_2000"] + 
    dataf$min_dis_cap750k_cntry * mdl_gain["dis_city_750k_near"] + 
    forest_cv * mdl_gain["per_forest_2000"] + 
    (forest_cv^2) * mdl_gain["per_forest_2000_2"] +
    (forest_cv^3) * mdl_gain["per_forest_2000_3"] + 
    (forest_cv^4) * mdl_gain["per_forest_2000_4"] +
    cont_fix + 
    bio_fix
  
  pred_gain <- ifelse(is.na(npv_refor_val$ind_var), 0, pred_gain)
  pred_gain <- exp(pred_gain)
  # change deforestation to zero if predicted deforestation is negative (does not occur)
  pred_gain <- ifelse(pred_gain < 0, 0, pred_gain)
  pred_gain <- ifelse(is.na(start_for_cv_bin), 0, pred_gain)
  # calculate forest cover
  pred_gain <- ifelse(pred_gain+forest_cv > 1, 1 - forest_cv, pred_gain)
  
  
  tau = 1
  tau_diff = 1
  counter = 0

  while(abs(tau_diff) > 0.005) {
    
    counter = counter + 1
    
    if(counter > 100){
      break
      print(paste0("No convergence in loop at", counter))
    }
    
    current_tau = (tau[counter]^0.195)
    
    # create a variable of character variables 
    bio_fix = paste0("F_biome_num", dataf$biome_num_f)
    bio_fix = mdl_loss[match(bio_fix, names(mdl_loss))]
    
    # create continent variable 
    cont_fix = paste0("F_cont_num_f", (ifelse(dataf$cont_var == "LatinAmerica", 1,
                                              ifelse(dataf$cont_var == "Africa", 2, 3))))
    cont_fix = mdl_loss[match(cont_fix, names(mdl_loss))]
    
    # calculate biome reforestation category 
    biome_npv_cat = biome_cat_df$biome_cat[match(dataf$biome_num, biome_cat_df$biome_num)]
    cont_cat = ifelse(dataf$country == "Brazil", 4,
                      ifelse(dataf$country == "Cambodia", 5,
                             ifelse(dataf$country == "Colombia", 6,
                                    ifelse(dataf$country == "Indonesia", 7,
                                           ifelse(dataf$country == "Liberia", 8,
                                                  ifelse(dataf$country == "Malaysia", 9,
                                                         ifelse(dataf$country == "Peru", 10,
                                                                ifelse(dataf$cont_var == "LatinAmerica", 1,
                                                                       ifelse(dataf$cont_var == "Africa", 2, 3)))))))))
    
    match_df = data.frame("uni_ID" = dataf$uni_ID,
                          "biome_npv_cat" = biome_npv_cat,
                          "cont_cat" = cont_cat, 
                          stringsAsFactors = FALSE)
    match_df = left_join(match_df, plant_per_df, by = c("cont_cat"= "contin_cat"))
    
    npv_refor_val = left_join(match_df, npv_carbon_refor_df, by = c("biome_npv_cat" = "biome_category",
                                                                    "cont_cat" = "contin_cat"))
    npv_refor_val$ind_var = ifelse(is.na(npv_refor_val$biome_npv_cat), 0, 1)
    
    
    bau_carbon_money = ((dataf[, carbon_name] * 3.67) + dataf$soil_crb + dataf$peat_crb) *
      (0 / 9.567)
    
    # calculate prediction model 
    bau_pred_loss = mdl_loss["(Intercept)"] + 
      mdl_loss[ag_rev_name] * ((dataf[,ag_rev_name] * 1) - bau_carbon_money)  +
      dataf$elev * mdl_loss["elev"] + 
      dataf$slope * mdl_loss["slope"] + 
      dataf$iucnlow_2000 * mdl_loss["iucnlow_2000"] + 
      dataf$iucnhigh_2000 * mdl_loss["iucnhigh_2000"] + 
      dataf$min_dis_cap750k_cntry * mdl_loss["dis_city_750k_near"] + 
      forest_cv * mdl_loss["per_forest_2000"] + 
      (forest_cv^2) * mdl_loss["per_forest_2000_2"] +
      (forest_cv^3) * mdl_loss["per_forest_2000_3"] + 
      (forest_cv^4) * mdl_loss["per_forest_2000_4"] +
      cont_fix + 
      bio_fix
    
    bau_pred_loss <- exp(bau_pred_loss)
    # change deforestation to zero if predicted deforestation is negative (does not occur)
    bau_pred_loss <- ifelse(bau_pred_loss < 0, 0, bau_pred_loss)
    bau_pred_loss <- ifelse(bau_pred_loss > forest_cv, forest_cv - 0, bau_pred_loss)
    bau_pred_loss_skm <- bau_pred_loss * cont_df$shp_size
    
    
    carbon_money = ((dataf[, carbon_name] * 3.67) + dataf$soil_crb + dataf$peat_crb) *
      (carbon_price / 9.567)
    
    
    # calculate prediction model 
    pred_loss = mdl_loss["(Intercept)"] + 
      mdl_loss[ag_rev_name] * ((dataf[,ag_rev_name] * current_tau) - carbon_money)  +
      dataf$elev * mdl_loss["elev"] + 
      dataf$slope * mdl_loss["slope"] + 
      dataf$iucnlow_2000 * mdl_loss["iucnlow_2000"] + 
      dataf$iucnhigh_2000 * mdl_loss["iucnhigh_2000"] + 
      dataf$min_dis_cap750k_cntry * mdl_loss["dis_city_750k_near"] + 
      forest_cv * mdl_loss["per_forest_2000"] + 
      (forest_cv^2) * mdl_loss["per_forest_2000_2"] +
      (forest_cv^3) * mdl_loss["per_forest_2000_3"] + 
      (forest_cv^4) * mdl_loss["per_forest_2000_4"] +
      cont_fix + 
      bio_fix
    
    pred_loss <- exp(pred_loss)
    # change deforestation to zero if predicted deforestation is negative (does not occur)
    pred_loss <- ifelse(pred_loss < 0, 0, pred_loss)
    pred_loss <- ifelse(pred_loss > forest_cv, forest_cv - 0, pred_loss)
    pred_loss_skm <- pred_loss * cont_df$shp_size

    def_bau =  sum(bau_pred_loss_skm, na.rm = TRUE) / sum(pred_loss_skm, na.rm = TRUE) 
    
    if(length(tau) > 1){
      
      tau[counter+1] = mean(c(tau[2:length(tau)], def_bau))
      
    } else{
      
      tau[counter+1] = def_bau
      
    }
    
    tau_diff = tau[counter+1] - tau[counter]
    
  }
  
  # calculate forest cover
  new_forest_cover = forest_cv - pred_loss + pred_gain
  
  return(list(new_forest_cover, pred_loss, pred_gain))  
  

}

############################################################################
# simulation function that estimates deforestation and emission
# with a carbon price
pred_fun_refor <- function(mdl_loss, 
                           mdl_gain,
                           dataf, 
                           forest_cv,
                           ag_rev_name, 
                           carbon_name, 
                           carbon_price,
                           biome_cat_df,
                           start_for_cv_bin,
                           npv_carbon_refor_df,
                           plant_per_df,
                           soil_carbon_refor_df,
                           npv_soil_carbon_refor_df){
  
  
  ########### loss prediction 
  # calculate carbon revenue 
  carbon_money = ((dataf[, carbon_name] * 3.67) + dataf$soil_crb + dataf$peat_crb) *
    (0 / 9.567)
  
  # create a variable of character variables 
  bio_fix = paste0("F_biome_num", dataf$biome_num_f)
  bio_fix = mdl_loss[match(bio_fix, names(mdl_loss))]
  
  # create continent variable 
  cont_fix = paste0("F_cont_num_f", (ifelse(dataf$cont_var == "LatinAmerica", 1,
                                            ifelse(dataf$cont_var == "Africa", 2, 3))))
  cont_fix = mdl_loss[match(cont_fix, names(mdl_loss))]
  
  # calculate biome reforestation category 
  biome_npv_cat = biome_cat_df$biome_cat[match(dataf$biome_num, biome_cat_df$biome_num)]
  cont_cat = ifelse(dataf$country == "Brazil", 4,
                    ifelse(dataf$country == "Cambodia", 5,
                           ifelse(dataf$country == "Colombia", 6,
                                  ifelse(dataf$country == "Indonesia", 7,
                                         ifelse(dataf$country == "Liberia", 8,
                                                ifelse(dataf$country == "Malaysia", 9,
                                                       ifelse(dataf$country == "Peru", 10,
                                                              ifelse(dataf$cont_var == "LatinAmerica", 1,
                                                                     ifelse(dataf$cont_var == "Africa", 2, 3)))))))))
  
  match_df = data.frame("uni_ID" = dataf$uni_ID,
                        "biome_npv_cat" = biome_npv_cat,
                        "cont_cat" = cont_cat, 
                        stringsAsFactors = FALSE)
  match_df = left_join(match_df, plant_per_df, by = c("cont_cat"= "contin_cat"))
  
  npv_refor_val = left_join(match_df, npv_carbon_refor_df, by = c("biome_npv_cat" = "biome_category",
                                                                  "cont_cat" = "contin_cat"))
  npv_refor_val$ind_var = ifelse(is.na(npv_refor_val$biome_npv_cat), 0, 1)
  
  
  # calculate prediction model 
  pred_loss = mdl_loss["(Intercept)"] + 
    mdl_loss[ag_rev_name] * (dataf[,ag_rev_name] - carbon_money)  +
    dataf$elev * mdl_loss["elev"] + 
    dataf$slope * mdl_loss["slope"] + 
    dataf$iucnlow_2000 * mdl_loss["iucnlow_2000"] + 
    dataf$iucnhigh_2000 * mdl_loss["iucnhigh_2000"] + 
    dataf$min_dis_cap750k_cntry * mdl_loss["dis_city_750k_near"] + 
    forest_cv * mdl_loss["per_forest_2000"] + 
    (forest_cv^2) * mdl_loss["per_forest_2000_2"] +
    (forest_cv^3) * mdl_loss["per_forest_2000_3"] + 
    (forest_cv^4) * mdl_loss["per_forest_2000_4"] +
    cont_fix + 
    bio_fix
  
  pred_loss <- exp(pred_loss)
  # change deforestation to zero if predicted deforestation is negative (does not occur)
  pred_loss <- ifelse(pred_loss < 0, 0, pred_loss)
  pred_loss <- ifelse(pred_loss > forest_cv, forest_cv - 0, pred_loss)
  
  
  #### gain prediction 
  # create a variable of character variables 
  bio_fix = paste0("F_biome_num", dataf$biome_num_f)
  bio_fix = mdl_gain[match(bio_fix, names(mdl_gain))]
  
  # create continent variable 
  cont_fix = paste0("F_cont_num_f", (ifelse(dataf$cont_var == "LatinAmerica", 1,
                                            ifelse(dataf$cont_var == "Africa", 2, 3))))
  cont_fix = mdl_gain[match(cont_fix, names(mdl_gain))]
  
  
  # create soil carbon reforestation data 
  soil_df = data.frame("uni_ID" = dataf$uni_ID,
                       "cont_cat" = cont_cat, 
                       stringsAsFactors = FALSE)
  
  soil_df = left_join(soil_df, npv_soil_carbon_refor_df, by = c("cont_cat" = "match"))
  
  # carbon gain 
  carbon_gain = ((npv_refor_val$discount_rate_10per * 1.26 * 3.67) + 
                   (dataf$soil_crb_ref * soil_df$discount_rate_10per)) * (carbon_price / 9.567)
  
  # calculate prediction model 
  pred_gain = mdl_gain["(Intercept)"] + 
    mdl_gain[ag_rev_name] * (dataf[,ag_rev_name] - carbon_gain)  +
    dataf$elev * mdl_gain["elev"] + 
    dataf$slope * mdl_gain["slope"] + 
    dataf$iucnlow_2000 * mdl_gain["iucnlow_2000"] + 
    dataf$iucnhigh_2000 * mdl_gain["iucnhigh_2000"] + 
    dataf$min_dis_cap750k_cntry * mdl_gain["dis_city_750k_near"] + 
    forest_cv * mdl_gain["per_forest_2000"] + 
    (forest_cv^2) * mdl_gain["per_forest_2000_2"] +
    (forest_cv^3) * mdl_gain["per_forest_2000_3"] + 
    (forest_cv^4) * mdl_gain["per_forest_2000_4"] +
    cont_fix + 
    bio_fix
  
  pred_gain <- ifelse(is.na(npv_refor_val$ind_var), 0, pred_gain)
  pred_gain <- exp(pred_gain)
  # change deforestation to zero if predicted deforestation is negative (does not occur)
  pred_gain <- ifelse(pred_gain < 0, 0, pred_gain)
  pred_gain <- ifelse(is.na(start_for_cv_bin), 0, pred_gain)
  # calculate forest cover
  pred_gain <- ifelse(pred_gain+forest_cv > 1, 1 - forest_cv, pred_gain)
  
  # calculate forest cover
  new_forest_cover = forest_cv - pred_loss + pred_gain
  
  return(list(new_forest_cover, pred_loss, pred_gain))
  
}

glob_model_pred <- function(mdl_loss, 
                            mdl_gain,
                            dataf, 
                            forest_cv,
                            year_steps,
                            ag_rev_name, 
                            carbon_name, 
                            carbon_price,
                            npv_carbon_refor_df,
                            biome_cat_df,
                            carbon_stock_dec_df,
                            plant_per_df,
                            soil_carbon_refor_df,
                            npv_soil_carbon_refor_df){
  

  
  #mdl_loss = cont_loss_list
  #mdl_gain = cont_gain_list
  #dataf = cont_df
  #forest_cv = start_for_cv
  #year_steps = seq(2020, 2050, 10)
  #ag_rev_name = "agrev_dec_max_00_10"
  #carbon_name = "carbon_baccini"
  #carbon_price = 20
  #npv_carbon_refor_df = npv_carbon_refor
  #biome_cat_df = biome_cat
  #carbon_stock_dec_df = carbon_stock_dec
  #plant_per_df = plant_per
  #soil_carbon_refor_df = soil_carbon_refor
  #npv_soil_carbon_refor_df = npv_soil_carbon_refor
 
  
  for_cv_list_defor = list()
  for_cv_list_refor = list()
  def_list = list()
  gain_list = list()
  
  for_2000 = dataf$per_all_forest + dataf$per_for_loss
  start_for_cv_bin = ifelse((start_for_cv == 0) & (for_2000 == 0), NA, 1)
  
  # loop through every decade
  for(yr in 1:length(year_steps)){
    
    
    if(yr == 1){
      
      # create prediction of forest gain and losss
      pred_list_defor = pred_fun_defor(mdl_loss = cont_loss_list, 
                           mdl_gain = cont_gain_list,
                           dataf = cont_df, 
                           forest_cv = start_for_cv,
                           ag_rev_name = ag_rev_name, 
                           carbon_name = carbon_name, 
                           carbon_price = 0,
                           start_for_cv_bin = start_for_cv_bin,
                           biome_cat_df = biome_cat_df,
                           npv_carbon_refor_df = npv_carbon_refor_df,
                           plant_per_df = plant_per_df,
                           soil_carbon_refor_df = soil_carbon_refor_df,
                           npv_soil_carbon_refor_df = npv_soil_carbon_refor_df)
      
      #### store in the list object
      for_cv_list_defor[[yr]] = pred_list_defor[[1]]
      def_list[[yr]] =  pred_list_defor[[2]]
      
      # create prediction of forest gain and losss
      pred_list_refor = pred_fun_refor(mdl_loss = cont_loss_list, 
                                          mdl_gain = cont_gain_list,
                                          dataf = cont_df, 
                                          forest_cv = start_for_cv,
                                          ag_rev_name = ag_rev_name, 
                                          carbon_name = carbon_name, 
                                          carbon_price = 0,
                                          biome_cat_df = biome_cat_df,
                                          start_for_cv_bin = start_for_cv_bin,
                                          npv_carbon_refor_df = npv_carbon_refor_df,
                                          plant_per_df = plant_per_df,
                                       soil_carbon_refor_df = soil_carbon_refor_df,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor_df)
      
      #### store in the list object
      for_cv_list_refor[[yr]] = pred_list_refor[[1]]
      gain_list[[yr]] = pred_list_refor[[3]]
      
     
    } else {
  
      # create prediction of forest gain and losss
      pred_list_defor = pred_fun_defor(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = for_cv_list_defor[[yr-1]],
                                       ag_rev_name = ag_rev_name, 
                                       carbon_name = carbon_name, 
                                       carbon_price = carbon_price,
                                       biome_cat_df = biome_cat_df,
                                       start_for_cv_bin = start_for_cv_bin,
                                       npv_carbon_refor_df = npv_carbon_refor_df,
                                       plant_per_df = plant_per_df,
                                       soil_carbon_refor_df = soil_carbon_refor_df,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor_df)
      
      #### store in the list object
      for_cv_list_defor[[yr]] = pred_list_defor[[1]]
      def_list[[yr]] =  pred_list_defor[[2]]
      
      # create prediction of forest gain and losss
      pred_list_refor = pred_fun_refor(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = for_cv_list_refor[[yr-1]],
                                       ag_rev_name = ag_rev_name, 
                                       carbon_name = carbon_name, 
                                       carbon_price = carbon_price,
                                       biome_cat_df = biome_cat_df,
                                       start_for_cv_bin = start_for_cv_bin,
                                       npv_carbon_refor_df = npv_carbon_refor_df,
                                       plant_per_df = plant_per_df,
                                       soil_carbon_refor_df = soil_carbon_refor_df,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor_df)
      
      #### store in the list object
      for_cv_list_refor[[yr]] = pred_list_refor[[1]]
      gain_list[[yr]] = pred_list_refor[[3]]
      
      
    }
    
  }
  
  ####################################################################
  # create data.frame
  # calculate biome reforestation category 
  biome_three_var = biome_cat_df$biome_cat[match(dataf$biome_num,
                                                 biome_cat_df$biome_num)]
  # create continent variable
  cont_cat = ifelse(dataf$cont_var == "LatinAmerica", 1,
                    ifelse(dataf$cont_var == "Africa", 2, 3))
  # create match data.frame 
  match_df = data.frame("uni_ID" = dataf$uni_ID,
                        "biome_car_cat" = biome_three_var,
                        "cont_cat" = cont_cat, 
                        stringsAsFactors = FALSE)
  # merge the dataset
  carb_stock_val = left_join(match_df, carbon_stock_dec_df, by = c("biome_car_cat" = "biome_category",
                                                                   "cont_cat" = "contin_cat"))
  carb_stock_val$ind_var = ifelse(is.na(carb_stock_val$biome_car_cat), 0, 1)
  
  
  # calculate biome reforestation category 
  biome_three_var = biome_cat_df$biome_cat[match(dataf$biome_num, biome_cat_df$biome_num)]
  cont_cat = ifelse(dataf$country == "Brazil", 4,
                    ifelse(dataf$country == "Cambodia", 5,
                           ifelse(dataf$country == "Colombia", 6,
                                  ifelse(dataf$country == "Indonesia", 7,
                                         ifelse(dataf$country == "Liberia", 8,
                                                ifelse(dataf$country == "Malaysia", 9,
                                                       ifelse(dataf$country == "Peru", 10,
                                                              ifelse(dataf$cont_var == "LatinAmerica", 1,
                                                                     ifelse(dataf$cont_var == "Africa", 2, 3)))))))))
  
  # create match data.frame 
  match_df = data.frame("uni_ID" = dataf$uni_ID,
                        "biome_car_cat" = biome_three_var,
                        "cont_cat" = cont_cat, 
                        stringsAsFactors = FALSE)
  match_df = left_join(match_df, plant_per_df, by = c("cont_cat"= "contin_cat"))
  
  carb_stock_val = left_join(match_df, carbon_stock_dec_df, by = c("biome_car_cat" = "biome_category",
                                                                   "cont_cat" = "contin_cat"))
  carb_stock_val$ind_var = ifelse(is.na(carb_stock_val$biome_car_cat), 0, 1)
  
  
  # create soil carbon reforestation data 
  soil_df = data.frame("uni_ID" = dataf$uni_ID,
                       "cont_cat" = cont_cat, 
                       stringsAsFactors = FALSE)
  
  soil_df = left_join(soil_df, soil_carbon_refor_df, by = c("cont_cat" = "match"))
  
  ####################################################################
  
  
  # calculate forest cover and deforestation in square kilometer
  # also calculate emissions for carbon price scenario and bau scenario
  for_cv_skm_defor_list <- lapply(for_cv_list_defor, function(x) x * cont_df$shp_size) 
  def_cv_skm_list <- lapply(def_list, function(x) x * cont_df$shp_size) 
  def_emiss_list <- lapply(def_list, function(x) x * cont_df$shp_size * ((cont_df[, carbon_name] * 3.67 * 100) + 
                                                                           (cont_df$soil_crb * 100) + (cont_df$peat_crb * 100)))  
  # gain square kilometer list
  for_cv_skm_refor_list <- lapply(for_cv_list_refor, function(x) x * cont_df$shp_size) 
  gain_cv_skm_list <- lapply(gain_list, function(x) x * cont_df$shp_size * carb_stock_val$ind_var) 
  gain_emiss_list <- list()
  
  for(gc in 1:length(year_steps)){
    
    carb_df = data.frame(carb_stock_val[, as.character(rev(1:gc))])
    gain_emiss_list[[gc]] <- apply(sapply(1:gc, function(x) gain_list[[x]] * 
                                            cont_df$shp_size * 
                                            100 * 
                                            ((carb_stock_val$ind_var *
                                                carb_df[,x] * 
                                                1.26 *
                                                3.67) + 
                                               (cont_df$soil_crb_ref * soil_df$discount_rate_10per))), 1, sum, na.rm = T) 
    
    ###+ (gain_list[[gc]] * cont_df$shp_size *  carb_stock_val$ind_var * 100 *
    ###     (dataf$soil_crb_ref * (match_df$nonplant * 0.21 + match_df$plant * 0.0085)))
    
    
  }
  
  
  # create a set of data frames 
  for_cv_skm_defor_df = as.data.frame(do.call(cbind, for_cv_skm_defor_list))
  def_cv_skm_df = as.data.frame(do.call(cbind, def_cv_skm_list))
  def_emiss_df = as.data.frame(do.call(cbind, def_emiss_list))
  for_cv_skm_refor_df = as.data.frame(do.call(cbind, for_cv_skm_refor_list))
  gain_cv_skm_df  = as.data.frame(do.call(cbind, gain_cv_skm_list))
  gain_emiss_df  = as.data.frame(do.call(cbind, gain_emiss_list))
  
  gain_df = (cont_df$per_for_gain * cont_df$shp_size * 100) * 
    ((1.26 *
        carb_stock_val$ind_var *
        3.67 * 
        carb_stock_val[, c("2", "3", "4", "5")]) +
       (cont_df$soil_crb_ref * soil_df$discount_rate_10per))
  colnames(gain_df) = colnames(gain_emiss_df)
  
  gain_emiss_df = data.frame(as.matrix(gain_emiss_df) + as.matrix(gain_df), stringsAsFactors = FALSE)
  
  
  # set column names 
  colnames(for_cv_skm_defor_df) = year_steps
  colnames(def_cv_skm_df) = year_steps
  colnames(def_emiss_df) = year_steps
  colnames(for_cv_skm_refor_df) = year_steps
  colnames(gain_cv_skm_df) = year_steps
  colnames(gain_emiss_df) = year_steps
  
  # create uni id and add continent and country 
  for_cv_skm_defor_df$uni_ID = cont_df$uni_ID
  for_cv_skm_defor_df$country = cont_df$country
  for_cv_skm_defor_df$continent = cont_df$cont_var
  
  def_cv_skm_df$uni_ID = cont_df$uni_ID
  def_cv_skm_df$country = cont_df$country
  def_cv_skm_df$continent = cont_df$cont_var
  
  def_emiss_df$uni_ID = cont_df$uni_ID
  def_emiss_df$country = cont_df$country
  def_emiss_df$continent = cont_df$cont_var
  
  for_cv_skm_refor_df$uni_ID = cont_df$uni_ID
  for_cv_skm_refor_df$country = cont_df$country
  for_cv_skm_refor_df$continent = cont_df$cont_var
  
  gain_cv_skm_df$uni_ID = cont_df$uni_ID
  gain_cv_skm_df$country = cont_df$country
  gain_cv_skm_df$continent = cont_df$cont_var
  
  gain_emiss_df$uni_ID = cont_df$uni_ID
  gain_emiss_df$country = cont_df$country
  gain_emiss_df$continent = cont_df$cont_var
  
  return(list(for_cv_skm_defor_df,
              def_cv_skm_df,
              def_emiss_df,
              for_cv_skm_refor_df,
              gain_cv_skm_df,
              gain_emiss_df))
  
  
  
  
}

library(dplyr)
library(tidyr)


###############################################################################
# set file folder
fold_path = "D:/Research/global_refor_mac/results/rds_reg_files/"

# upload simulation regressions and dataframe for all three continents
cont_loss_list <- readRDS(paste0(fold_path, "reg_loss_cont_v1.rds"))
cont_loss_list <- cont_loss_list$coefficients
cont_loss_list[is.na(cont_loss_list)] <- 0
# upload simulation regressions and dataframe for all three continents
cont_gain_list <- readRDS(paste0(fold_path, "reg_gain_cont_v1.rds"))
cont_gain_list <- cont_gain_list$coefficients
cont_gain_list[is.na(cont_gain_list)] <- 0
cont_df <-  readRDS(paste0(fold_path, "all_cont_reg.rds"))
cont_df$carbon_baccini <- cont_df$carbon_baccini * 1.26 / 2 # adjust to carbon and divide by two 
cont_df$peat_crb <- (cont_df$peat_crb / 815.545) * 1782 # change peat CO2 to "new" peat value
cont_df$soil_crb_ref <- (1 / 0.085) * cont_df$soil_crb

# read in biome category dataset
biome_cat = read.csv("D:/Research/global_refor_mac/data/carbon_file/biome_cat.csv",
                     stringsAsFactors = FALSE)
npv_carbon_refor = read.csv("D:/Research/global_refor_mac/data/carbon_file/npv_carbon_refor_trans_new.csv",
                            stringsAsFactors = FALSE)

# read in the carbon stock per decade
carbon_stock_dec = read.csv("D:/Research/global_refor_mac/data/carbon_file/carbon_stock_per_decade_new.csv",
                            stringsAsFactors = FALSE)
colnames(carbon_stock_dec)[4:13] = 1:10

# read in the carbon stock per decade
plant_per = read.csv("D:/Research/global_refor_mac/data/carbon_file/plant_perc.csv",
                     stringsAsFactors = FALSE)


# read in soil npv carbon refor
soil_carbon_refor = read.csv("D:/Research/global_refor_mac/data/carbon_file/soil_carbon_stock_per_decade.csv",
                             stringsAsFactors = FALSE)
npv_soil_carbon_refor = read.csv("D:/Research/global_refor_mac/data/carbon_file/npv_soil_carbon_stock_per_decade.csv",
                                 stringsAsFactors = FALSE)



##################################################################################
#### create bau scenario 

# adjust the data.frame
cont_df = cont_df[which(is.na(cont_df$carbon_baccini) == FALSE), ]
cont_df$soil_crb[is.na(cont_df$soil_crb)] <- 0
cont_df$peat_crb[is.na(cont_df$peat_crb)] <- 0
start_for_cv = cont_df$per_all_forest + cont_df$per_for_gain


#### Zero dollars
glob_model_0dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                      mdl_gain = cont_gain_list,
                                      dataf = cont_df, 
                                      forest_cv = start_for_cv,
                                      year_steps = seq(2020, 2050, 10),
                                      ag_rev_name = "agrev_dec_max_00_10", 
                                      carbon_name = "carbon_baccini", 
                                      carbon_price = 0,
                                      npv_carbon_refor_df = npv_carbon_refor,
                                      biome_cat_df = biome_cat,
                                      carbon_stock_dec_df = carbon_stock_dec,
                                      plant_per_df = plant_per,
                                      soil_carbon_refor_df = soil_carbon_refor,
                                      npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_5dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                      mdl_gain = cont_gain_list,
                                      dataf = cont_df, 
                                      forest_cv = start_for_cv,
                                      year_steps = seq(2020, 2050, 10),
                                      ag_rev_name = "agrev_dec_max_00_10", 
                                      carbon_name = "carbon_baccini", 
                                      carbon_price = 5,
                                      npv_carbon_refor_df = npv_carbon_refor,
                                      biome_cat_df = biome_cat,
                                      carbon_stock_dec_df = carbon_stock_dec,
                                      plant_per_df = plant_per,
                                      soil_carbon_refor_df = soil_carbon_refor,
                                      npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_10dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 10,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_20dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 20,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_30dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 30,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_40dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 40,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_50dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 50,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_60dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 60,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)


#### Zero dollars
glob_model_70dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 70,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_80dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 80,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

#### Zero dollars
glob_model_90dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                       mdl_gain = cont_gain_list,
                                       dataf = cont_df, 
                                       forest_cv = start_for_cv,
                                       year_steps = seq(2020, 2050, 10),
                                       ag_rev_name = "agrev_dec_max_00_10", 
                                       carbon_name = "carbon_baccini", 
                                       carbon_price = 90,
                                       npv_carbon_refor_df = npv_carbon_refor,
                                       biome_cat_df = biome_cat,
                                       carbon_stock_dec_df = carbon_stock_dec,
                                       plant_per_df = plant_per,
                                       soil_carbon_refor_df = soil_carbon_refor,
                                       npv_soil_carbon_refor_df = npv_soil_carbon_refor)

glob_model_100dollar <- glob_model_pred(mdl_loss = cont_loss_list, 
                                        mdl_gain = cont_gain_list,
                                        dataf = cont_df, 
                                        forest_cv = start_for_cv,
                                        year_steps = seq(2020, 2050, 10),
                                        ag_rev_name = "agrev_dec_max_00_10", 
                                        carbon_name = "carbon_baccini", 
                                        carbon_price = 100,
                                        npv_carbon_refor_df = npv_carbon_refor,
                                        biome_cat_df = biome_cat,
                                        carbon_stock_dec_df = carbon_stock_dec,
                                        plant_per_df = plant_per,
                                        soil_carbon_refor_df = soil_carbon_refor,
                                        npv_soil_carbon_refor_df = npv_soil_carbon_refor)


############################################################################
### save the model out 
out_path = "D:/Research/global_refor_mac/sensitivity_runs/"

#### create list object data.frame
model_list = list(glob_model_0dollar, glob_model_5dollar, 
                  glob_model_10dollar, glob_model_20dollar, 
                  glob_model_30dollar, glob_model_40dollar, 
                  glob_model_50dollar, glob_model_60dollar, 
                  glob_model_70dollar, glob_model_80dollar, 
                  glob_model_90dollar, glob_model_100dollar)

saveRDS(model_list, paste0(out_path, "refor_leakage_model_optionb.rds"))

#for_cv_skm_df,
#def_cv_skm_df,
#def_emiss_df,
#gain_cv_skm_df,
#gain_emiss_df

